unit Projection; {************************************************************************} {* Projection.p *} {* by Michael Castle (Pascal) and Janice Keller (assembly code) *} {* University of Michigan Mental Health Research Institute (MHRI) *} {* e-mail address: mike.castle@med.umich.edu *} {* ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * } interface uses Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils,Resources, Errors, Palettes, Globals, Utilities, File2, File1, Graphics, Camera, Filters, Stacks; procedure DoProjection; function ShowProjectDialogBox: boolean; procedure Project; implementation const bigpowerof2 = 8192; {used for integer trigonometric arithmetic} type MHRIRealArray = array[1..MaxPics] of extended; RealPoint = record x: extended; y: extended; end; {record} IntPtr = ^integer; LongPtr = ^LongInt; var SliceInterval: extended; {distance, in pixels, between original slices} BoundRect: rect; {boundary rectangle for a rectangular selection} cancelled: boolean; ProjSize: LongInt; {******************************************************************************} {* This procedure frees memory allocated for buffers used in projection calculations. *} {******************************************************************************} procedure DisposeProjectionPtrs (projptr, opaptr, brightcueptr: ptr; zbufptr, countbufptr, cuezbufptr: IntPtr; sumbufptr: LongPtr); begin if zbufptr <> nil then DisposePtr(Ptr(zbufptr)); if opaptr <> nil then DisposePtr(opaptr); if sumbufptr <> nil then DisposePtr(Ptr(sumbufptr)); if countbufptr <> nil then DisposePtr(Ptr(countbufptr)); if cuezbufptr <> nil then DisposePtr(Ptr(cuezbufptr)); if brightcueptr <> nil then DisposePtr(brightcueptr); if projptr <> nil then DisposePtr(projptr); end; {******************************************************************************} {* This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about *} {* the x-axis. Integer arithmetic, precomputation of values, and iterative addition rather than multiplication *} {* inside a loop are used extensively to make the code run efficiently. Projection parameters stored in global *} {* variables determine how the projection will be performed. This procedure returns various buffers which *} {* are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*} {******************************************************************************} procedure DoOneProjectionX (nSlices, ycenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string); var i, j, k, {loop indices} thispixel: integer; {the current pixel to be projected} offset, offsetinit: longint; {precomputed offsets into an image buffer} projaddr, {buffer address for final projected image} opaaddr, {buffer address for opacity surface projection component} brightcueaddr, {buffer address for brightest-point projection with interior depth cues} zbufaddr, {z-buffer address for surface projections (nearest-point/opacity)} cuezbufaddr, {z-buffer address for brightest-point projection with interior depth cues} countbufaddr, {buffer addr for counting pixels in mean-value projection} sumbufaddr: longint; {buffer addr for summing pixels in mean-value projection} z, {z-coordinate of points in current slice before rotation} ynew, znew, {y- and z-coordinates of current point after rotation} zmax, zmin, {z-coordinates of first and last slices before rotation} zmaxminuszmintimes100: longint; {precomputed values to save time in loops} c100minusDepthCueInt, c100minusDepthCueSurf: integer; DepthCueIntlessthan100, DepthCueSurflessthan100: boolean; OpacityOrNearestPt, OpacityAndNotNearestPt: boolean; MeanVal, BrightestPt: boolean; ysintheta, ycostheta, {values used in rotational calculations before projection} zsintheta, zcostheta, ysinthetainit, ycosthetainit: longint; p, {pointer to final projected image buffer} op, {pointer to opacity surface projection buffer} bp: ptr; {pointer to brightest-point projection buffer with interior depth cues} zp, {pointer to surface projection (nearest-point/opacity) z-buffer} qp, {pointer to z-buffer for brightest-point projection with interior depth cues} cp: IntPtr; {pointer to buffer for counting pixels for mean-value projection} sp: LongPtr; {pointer to mean-value summing buffer} width: integer; theLine: LineType; begin with BoundRect do begin {precompute values to avoid computations in loop} width := right - left; zmax := zcenter + projheight div 2; {find z-coordinates of first and last slices} zmin := zcenter - projheight div 2; zmaxminuszmintimes100 := 100 * (zmax - zmin); c100minusDepthCueInt := 100 - DepthCueInt; c100minusDepthCueSurf := 100 - DepthCueSurf; DepthCueIntlessthan100 := DepthCueInt < 100; DepthCueSurflessthan100 := DepthCueSurf < 100; OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0); OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint); MeanVal := ProjectionMethod = MeanValue; BrightestPt := ProjectionMethod = BrightestPoint; projaddr := ord4(projptr); opaaddr := ord4(opaptr); brightcueaddr := ord4(brightcueptr); zbufaddr := ord4(zbufptr); cuezbufaddr := ord4(cuezbufptr); countbufaddr := ord4(countbufptr); sumbufaddr := ord4(sumbufptr); ycosthetainit := (top - ycenter - 1) * costheta; ysinthetainit := (top - ycenter - 1) * sintheta; offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 - 1; for k := 1 to nSlices do begin SelectSlice(k); z := round((k - 1) * SliceInterval) - zcenter; zcostheta := z * costheta; zsintheta := z * sintheta; ycostheta := ycosthetainit; ysintheta := ysinthetainit; for j := top to bottom - 1 do begin {read each row in the current slice} ycostheta := ycostheta + costheta; {rotate about x-axis and find new y,z} ysintheta := ysintheta + sintheta; {x-coordinates will not change} ynew := (ycostheta - zsintheta) div bigpowerof2 + ycenter - top; znew := (ysintheta + zcostheta) div bigpowerof2 + zcenter; offset := offsetinit + ynew * projwidth; GetLine(left, j, width, theLine); for i := 0 to width - 1 do begin {read each pixel in current row and project it} thispixel := theLine[i]; offset := offset + 1; if (offset >= ProjSize) or (offset < 0) then offset := 0; if (thispixel <= TransparencyUpper) and (thispixel >= TransparencyLower) then begin if OpacityOrNearestPt then begin zp := IntPtr(zbufaddr + offset + offset); if (znew < zp^) then begin zp^ := znew; if OpacityAndNotNearestPt then begin op := ptr(opaaddr + offset); if (DepthCueSurflessthan100) then op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100) else op^ := thispixel; end else begin p := ptr(projaddr + offset); if DepthCueSurflessthan100 then p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100) else p^ := thispixel; end; end; end; {NearestPoint case} if MeanVal then begin sp := LongPtr(sumbufaddr + offset + offset + offset + offset); sp^ := sp^ + thispixel; cp := IntPtr(countbufaddr + offset + offset); cp^ := cp^ + 1; end {MeanValue case} else if BrightestPt then begin if (DepthCueIntlessthan100) then begin p := ptr(projaddr + offset); bp := ptr(brightcueaddr + offset); qp := IntPtr(cuezbufaddr + offset + offset); if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (znew < qp^)) then begin bp^ := thispixel; {use z-buffer to ensure that if depth-cueing is on, } qp^ := znew; {the closer of two equally-bright points is displayed} p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100); end; {then} end else begin p := ptr(projaddr + offset); if (thispixel < BAND(p^, 255)) then p^ := thispixel; end; {else} end; {BrightestPoint case} end; end; {for j} end; {for i} UpdateMeter(10 + (k * 90) div nSlices, str); if CommandPeriod then begin cancelled := true; beep; leave; end; end; {for k} end; {with} end; {procedure DoOneProjectionX} {******************************************************************************} {* This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about *} {* the y-axis. Integer arithmetic, precomputation of values, and iterative addition rather than multiplication *} {* inside a loop are used extensively to make the code run efficiently. Projection parameters stored in global *} {* variables determine how the projection will be performed. This procedure returns various buffers which *} {* are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*} {******************************************************************************} procedure DoOneProjectionY (nSlices, xcenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string); var i, j, k, thispixel: integer; offset, offsetinit: longint; projaddr, opaaddr, brightcueaddr, zbufaddr, cuezbufaddr, countbufaddr, sumbufaddr: longint; z, xnew, znew, zmax, zmin, zmaxminuszmintimes100: longint; c100minusDepthCueInt, c100minusDepthCueSurf: integer; DepthCueIntlessthan100, DepthCueSurflessthan100: boolean; OpacityOrNearestPt, OpacityAndNotNearestPt: boolean; MeanVal, BrightestPt: boolean; xsintheta, xcostheta, zsintheta, zcostheta, xsinthetainit, xcosthetainit: longint; p, op, bp: ptr; zp, qp, cp: IntPtr; sp: LongPtr; width: integer; aLine: LineType; begin with BoundRect do begin width := right - left; zmax := zcenter + projwidth div 2; zmin := zcenter - projwidth div 2; zmaxminuszmintimes100 := 100 * (zmax - zmin); c100minusDepthCueInt := 100 - DepthCueInt; c100minusDepthCueSurf := 100 - DepthCueSurf; DepthCueIntlessthan100 := DepthCueInt < 100; DepthCueSurflessthan100 := DepthCueSurf < 100; OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0); OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint); MeanVal := ProjectionMethod = MeanValue; BrightestPt := ProjectionMethod = BrightestPoint; projaddr := ord4(projptr); opaaddr := ord4(opaptr); brightcueaddr := ord4(brightcueptr); zbufaddr := ord4(zbufptr); cuezbufaddr := ord4(cuezbufptr); countbufaddr := ord4(countbufptr); sumbufaddr := ord4(sumbufptr); xcosthetainit := (left - xcenter - 1) * costheta; xsinthetainit := (left - xcenter - 1) * sintheta; for k := 1 to nSlices do begin SelectSlice(k); z := round((k - 1) * SliceInterval) - zcenter; zcostheta := z * costheta; zsintheta := z * sintheta; offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 - projwidth; for j := top to bottom - 1 do begin xcostheta := xcosthetainit; xsintheta := xsinthetainit; offsetinit := offsetinit + projwidth; GetLine(left, j, width, aLine); for i := 0 to width - 1 do begin thispixel := aLine[i]; xcostheta := xcostheta + costheta; xsintheta := xsintheta + sintheta; if (thispixel <= TransparencyUpper) and (thispixel >= TransparencyLower) then begin xnew := (xcostheta + zsintheta) div bigpowerof2 + xcenter - left; znew := (zcostheta - xsintheta) div bigpowerof2 + zcenter; offset := offsetinit + xnew; if (offset >= ProjSize) or (offset < 0) then offset := 0; if OpacityOrNearestPt then begin zp := IntPtr(zbufaddr + offset + offset); if (znew < zp^) then begin zp^ := znew; if OpacityAndNotNearestPt then begin op := ptr(opaaddr + offset); if (DepthCueSurflessthan100) then op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100) else op^ := thispixel; end else begin p := ptr(projaddr + offset); if DepthCueSurflessthan100 then p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100) else p^ := thispixel; end; end; end; {NearestPoint case} if MeanVal then begin sp := LongPtr(sumbufaddr + offset + offset + offset + offset); sp^ := sp^ + thispixel; cp := IntPtr(countbufaddr + offset + offset); cp^ := cp^ + 1; end {MeanValue case} else if BrightestPt then begin if (DepthCueIntlessthan100) then begin p := ptr(projaddr + offset); bp := ptr(brightcueaddr + offset); qp := IntPtr(cuezbufaddr + offset + offset); if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (znew < qp^)) then begin bp^ := thispixel; qp^ := znew; p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100); end; {then} end else begin p := ptr(projaddr + offset); if (thispixel < BAND(p^, 255)) then p^ := thispixel; end; {else} end; {BrightestPoint case} end; end; {for j} end; {for i} UpdateMeter(10 + (k * 90) div nSlices, str); if CommandPeriod then begin cancelled := true; beep; leave; end; end; {for k} end; {with} end; {procedure DoOneProjectionY} {******************************************************************************} {* This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about *} {* the z-axis. Integer arithmetic, precomputation of values, and iterative addition rather than multiplication *} {* inside a loop are used extensively to make the code run efficiently. Projection parameters stored in global *} {* variables determine how the projection will be performed. This procedure returns various buffers which *} {* are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*} {******************************************************************************} procedure DoOneProjectionZ (nSlices, xcenter, ycenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string); var i, j, k, thispixel: integer; offset, offsetinit: longint; projaddr, opaaddr, brightcueaddr, zbufaddr, cuezbufaddr, countbufaddr, sumbufaddr: longint; z, xnew, ynew, zmax, zmin, zmaxminuszmintimes100: longint; c100minusDepthCueInt, c100minusDepthCueSurf: integer; DepthCueIntlessthan100, DepthCueSurflessthan100: boolean; OpacityOrNearestPt, OpacityAndNotNearestPt: boolean; MeanVal, BrightestPt: boolean; xsintheta, xcostheta, ysintheta, ycostheta: longint; xsinthetainit, xcosthetainit, ysinthetainit, ycosthetainit: longint; p, op, bp: ptr; zp, qp, cp: IntPtr; sp: LongPtr; width: integer; theLine: LineType; begin with BoundRect do begin width := right - left; zmax := zcenter + projwidth div 2; zmin := zcenter - projwidth div 2; zmaxminuszmintimes100 := 100 * (zmax - zmin); c100minusDepthCueInt := 100 - DepthCueInt; c100minusDepthCueSurf := 100 - DepthCueSurf; DepthCueIntlessthan100 := DepthCueInt < 100; DepthCueSurflessthan100 := DepthCueSurf < 100; OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0); OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint); MeanVal := ProjectionMethod = MeanValue; BrightestPt := ProjectionMethod = BrightestPoint; projaddr := ord4(projptr); opaaddr := ord4(opaptr); brightcueaddr := ord4(brightcueptr); zbufaddr := ord4(zbufptr); cuezbufaddr := ord4(cuezbufptr); countbufaddr := ord4(countbufptr); sumbufaddr := ord4(sumbufptr); xcosthetainit := (left - xcenter - 1) * costheta; xsinthetainit := (left - xcenter - 1) * sintheta; ycosthetainit := (top - ycenter - 1) * costheta; ysinthetainit := (top - ycenter - 1) * sintheta; offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 + left - 1; for k := 1 to nSlices do begin SelectSlice(k); z := round((k - 1) * SliceInterval) - zcenter; ycostheta := ycosthetainit; ysintheta := ysinthetainit; for j := top to bottom - 1 do begin ycostheta := ycostheta + costheta; ysintheta := ysintheta + sintheta; xcostheta := xcosthetainit; xsintheta := xsinthetainit; GetLine(left, j, width, theLine); for i := 0 to width - 1 do begin thisPixel := theLine[i]; xcostheta := xcostheta + costheta; xsintheta := xsintheta + sintheta; if (thispixel <= TransparencyUpper) and (thispixel >= TransparencyLower) then begin xnew := (xcostheta - ysintheta) div bigpowerof2 + xcenter - left; ynew := (xsintheta + ycostheta) div bigpowerof2 + ycenter - top; offset := offsetinit + ynew * projwidth + xnew; if (offset >= ProjSize) or (offset < 0) then offset := 0; if OpacityOrNearestPt then begin zp := IntPtr(zbufaddr + offset + offset); if (z < zp^) then begin zp^ := z; if OpacityAndNotNearestPt then begin op := ptr(opaaddr + offset); if (DepthCueSurflessthan100) then op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100) else op^ := thispixel; end else begin p := ptr(projaddr + offset); if DepthCueSurflessthan100 then p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100) else p^ := thispixel; end; end; end; {NearestPoint case} if MeanVal then begin sp := LongPtr(sumbufaddr + offset + offset + offset + offset); sp^ := sp^ + thispixel; cp := IntPtr(countbufaddr + offset + offset); cp^ := cp^ + 1; end {MeanValue case} else if BrightestPt then begin if (DepthCueIntlessthan100) then begin p := ptr(projaddr + offset); bp := ptr(brightcueaddr + offset); qp := IntPtr(cuezbufaddr + offset + offset); if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (z < qp^)) then begin bp^ := thispixel; qp^ := z; p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100); end; {then} end else begin p := ptr(projaddr + offset); if (thispixel < BAND(p^, 255)) then p^ := thispixel; end; {else} end; {BrightestPoint case} end; end; {for j} end; {for i} UpdateMeter(10 + (k * 90) div nSlices, str); if CommandPeriod then begin cancelled := true; beep; leave; end; end; {for k} end; {with} end; {procedure DoOneProjectionZ} {******************************************************************************} {* This code initializes buffers by filling them with uniform values. The *} {* The buffer can be filled with one, two or four byte values. *} {******************************************************************************} procedure InitializeBuffer (p: ptr; size: longint; value, step: integer); type IntPtrType=^integer; LongPtrType=^LongInt; var i, Lstep: longint; IntPtr:IntPtrTYpe; LongPtr:LongPtrType; begin Lstep:=step; case step of 1: for i := 1 to size do begin p^ := value; p := Ptr(ord4(p) + Lstep); end; 2: begin IntPtr:=IntPtrType(p); for i := 1 to size do begin IntPtr^ := value; IntPtr := IntPtrType(ord4(IntPtr) + Lstep); end; end; 4: begin LongPtr:=LongPtrType(p); for i := 1 to size do begin LongPtr^ := value; LongPtr := LongPtrType(ord4(LongPtr) + Lstep); end; end; end; end; {******************************************************************************} {* This procedure creates a sequence of projections of a rotating volume (stack of slices) onto a plane using *} {* nearest-point (surface), brightest-point, or mean-value projection or a weighted combination of nearest- *} {* point projection with either of the other two methods (partial opacity). The user may choose to rotate the *} {* volume about any of the three orthogonal axes (x,y, or z), make portions of the volume transparent, or add *} {* a greater degree of visual realism by employing depth cues. *} {******************************************************************************} procedure DoProjection; var tport: Grafptr; nSlices: integer; {number of slices in volume} projwidth, projheight: LongInt; {dimensions of projection image} xcenter, ycenter, zcenter: integer; {coordinates of center of volume of rotation} theta: integer; {current angle of rotation in degrees} thetarad: extended; {current angle of rotation in radians} sintheta, costheta: longint; {sine and cosine of current angle} xsinthetainit, xcosthetainit: longint; {precomputed products to avoid mult in loops} offset, MemoryRequired: longint; p, op, bp, projptr, opaptr, brightcueptr: ptr; zp, zbufptr, qp, cuezbufptr, countbufptr, cp: IntPtr; sp, sumbufptr: LongPtr; curval, prevval, nextval, aboveval, belowval: integer; ignore: integer; {irrelevant return value from a function} ticksinit, tickstogo, TicksForOneProjection: longint; str, TimeStr, seconds: str255; SaveProjectionsTemp, AutoSelectAll, AllocatingBuffers: boolean; n, nProjections, angle: integer; SourceInfo, DestInfo: InfoPtr; width, height: LongInt; procedure Abort; begin DisposeProjectionPtrs(projptr, opaptr, brightcueptr, zbufptr, countbufptr, cuezbufptr, sumbufptr); if AllocatingBuffers and (MaxBlock > 20000) then PutError('Insufficient Memory.'); AbortMacro; {exit(DoProjection);} {ppc-bug} end; begin ShowWatch; AutoSelectAll := not Info^.RoiShowing; if AutoSelectAll then SelectAll(false); if NotInBounds then exit(DoProjection); cancelled := false; SourceInfo := Info; GetPort(tPort); with Info^ do begin SetPort(GrafPtr(osPort)); BoundRect := Roirect; end; if (AngleInc = 0) and (TotalAngle <> 0) then AngleInc := 5; angle := 0; nProjections := 0; if AngleInc = 0 then nProjections := 1 else while angle <= TotalAngle do begin nProjections := nProjections + 1; angle := angle + AngleInc; end; if angle > 360 then nProjections := nProjections - 1; if nProjections <= 0 then nProjections := 1; nSlices := Info^.StackInfo^.nSlices; {get number of slices in volume} SliceInterval := info^.StackInfo^.SliceSpacing; with BoundRect do begin width := right - left; height := bottom - top; xcenter := (left + right) div 2; {find center of volume of rotation} ycenter := (top + bottom) div 2; zcenter := round(nSlices * SliceInterval / 2.0); if MinProjSize and (AxisOfRotation <> ZAxis) then begin case AxisOfRotation of {find dimensions of projection image} XAxis: begin projheight := round(sqrt(sqr(nSlices * SliceInterval) + height * height)); projwidth := right - left; end; {XAxis} YAxis: begin projwidth := round(sqrt(sqr(nSlices * SliceInterval) + ord4(right - left) * (right - left))); projheight := bottom - top; end; {YAxis} end; {case} end {then} else begin projwidth := round(sqrt(sqr(nSlices * SliceInterval) + width * width)); projheight := round(sqrt(sqr(nSlices * SliceInterval) + height * height)); end; {else make all windows the same size regardless of rotation axis} end; {with BoundRect} if odd(projwidth) then projwidth := projwidth + 1; projptr := nil; zbufptr := nil; opaptr := nil; brightcueptr := nil; cuezbufptr := nil; sumbufptr := nil; countbufptr := nil; AllocatingBuffers := true; projsize := projwidth * projheight; projptr := NewPtr(projsize); if projptr = nil then begin Abort; exit(DoProjection) end; if (ProjectionMethod = NearestPoint) or (Opacity > 0) then begin {get other required buffers} zbufptr := IntPtr(NewPtr(projsize * 2)); if zbufptr = nil then begin Abort; exit(DoProjection) end; end; if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin opaptr := NewPtr(projsize); if opaptr = nil then begin Abort; exit(DoProjection) end; end; if (ProjectionMethod = BrightestPoint) and (DepthCueInt < 100) then begin brightcueptr := NewPtr(projsize); cuezbufptr := IntPtr(NewPtr(projsize * 2)); if (brightcueptr = nil) or (cuezbufptr = nil) then begin Abort; exit(DoProjection) end; end; if (ProjectionMethod = MeanValue) then begin sumbufptr := LongPtr(NewPtr(projsize * 4)); countbufptr := IntPtr(NewPtr(projsize * 2)); if (sumbufptr = nil) or (countbufptr = nil) then begin Abort; exit(DoProjection) end; end; AllocatingBuffers := false; SaveProjectionsTemp := FALSE; {check whether we have enough memory to open} MemoryRequired := nProjections * projsize + SizeOf(PicInfo) + MinFree; if (MemoryRequired > FreeMem) and not (SaveProjections) then begin str := 'Insufficient memory to create entire stack of projections. Projection images will be saved to disk.'; if (PutMessageWithCancel(str) = cancel) then begin Abort; exit(DoProjection) end; SaveProjections := TRUE; SaveProjectionsTemp := TRUE; end; if (SaveProjections) then begin {prepare to save projections as created if desired} SaveAsWhat := AsTiff; SaveAllState := SaveAllStage1; end; TimeStr := '?'; theta := InitAngle; {begin rotation with user-specified angle} ticksinit := TickCount; for n := 1 to nProjections do begin if (SaveProjections) or (n = 1) then begin {open new window or stack slice} if SaveProjections then case AxisOfRotation of XAxis: str := StringOf('Projection X ', theta:3); YAxis: str := StringOf('Projection Y ', theta:3); ZAxis: str := StringOf('Projection Z ', theta:3); end else str := 'Projections'; if not NewPicWindow(str, projwidth, projheight) then begin Abort; exit(DoProjection) end; end else if not (AddSlice(false)) then begin Abort; exit(DoProjection) end; str := concat('Projecting: ', long2str(n), '/', long2str(nProjections), ' (', long2str(Theta), '¡', ', ', TimeStr, ')'); ShowMeter; UpdateMeter(0, str); thetarad := theta * pi / 180.0; costheta := round(bigpowerof2 * cos(thetarad)); sintheta := round(bigpowerof2 * sin(thetarad)); p := projptr; {initialize required projection buffers} InitializeBuffer(p, projsize, 255, 1); if (ProjectionMethod = NearestPoint) or (Opacity > 0) then begin zp := zbufptr; InitializeBuffer(Ptr(zp), projsize, 32767, 2); end; {then} if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin op := opaptr; InitializeBuffer(op, projsize, 255, 1); end; {then} if (ProjectionMethod = MeanValue) then begin sp := sumbufptr; cp := countbufptr; InitializeBuffer(Ptr(sp), projsize, 0, 4); InitializeBuffer(Ptr(cp), projsize, 0, 2); end; if (ProjectionMethod = BrightestPoint) and (DepthCueInt < 100) then begin bp := brightcueptr; InitializeBuffer(bp, projsize, 255, 1); qp := cuezbufptr; InitializeBuffer(Ptr(qp), projsize, 255, 2); end; {then} UpdateMeter(10, str); DestInfo := Info; Info := SourceInfo; case AxisOfRotation of XAxis: DoOneProjectionX(nSlices, ycenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str); YAxis: DoOneProjectionY(nSlices, xcenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str); ZAxis: DoOneProjectionZ(nSlices, xcenter, ycenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str); end; Info := DestInfo; if n = 1 then TicksForOneProjection := TickCount - TicksInit; TicksToGo := (nProjections - n) * TicksForOneProjection; NumToString((TicksToGo div 60) mod 60, seconds); if length(seconds) = 1 then seconds := concat('0', seconds); timestr := concat(long2str(TicksToGo div 3600), ':', seconds); if (ProjectionMethod = MeanValue) then begin p := projptr; {calculate the mean-value image from returned info} sp := sumbufptr; cp := countbufptr; for offset := 0 to projsize - 1 do begin if (cp^ > 0) then p^ := sp^ div cp^; p := ptr(ord4(p) + 1); sp := LongPtr(ord4(sp) + 4); cp := IntPtr(ord4(cp) + 2); end; {for} end; {then} if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin p := projptr; {calculate surface proj component (opacity on)} op := opaptr; {and combine with another proj component} for offset := 0 to projsize - 1 do begin p^ := (Opacity * BAND(op^, 255) + (100 - Opacity) * BAND(p^, 255)) div 100; p := ptr(ord4(p) + 1); op := ptr(ord4(op) + 1); end; {for} end; {then} if (AxisOfRotation = ZAxis) then begin {interpolate for z-axis rotation} p := ptr(ord4(projptr) + projwidth); for offset := projwidth to projsize - 1 - projwidth do begin curval := BAND(p^, 255); prevval := BAND(ptr(ord4(p) - 1)^, 255); nextval := BAND(ptr(ord4(p) + 1)^, 255); aboveval := BAND(ptr(ord4(p) - projwidth)^, 255); belowval := BAND(ptr(ord4(p) + projwidth)^, 255); if (curval = 255) and (prevval <> 255) and (nextval <> 255) and (aboveval <> 255) and (belowval <> 255) then p^ := (prevval + nextval + aboveval + belowval) div 4; p := ptr(ord4(p) + 1); end; end; if (SaveProjections) or (n = 1) then BlockMove(projptr, Info^.PicBaseAddr, projsize) {whole ROI write to projection image} else BlockMove(projptr, Info^.StackInfo^.PicBaseH[n]^, projsize); UpdateMeter(-1, ''); {dispose of meter} if cancelled then begin if n > 1 then DeleteSlice; leave; end; if (SaveProjections) then begin SaveAs('', 0); {just save and put away current image after creating it} ignore := CloseAWindow(info^.wptr); end else if n = 1 then begin {create new stack for first projection if not saving projections} if not MakeStackFromWindow then begin Abort; exit(DoProjection) end end; theta := (theta + AngleInc) mod 360; UpdatePicWindow; end; {for} SaveAllState := NoSaveAll; if SaveProjectionsTemp then {turn this back off if we turned it on due to lack of memory} SaveProjections := FALSE; DisposeProjectionPtrs(projptr, opaptr, brightcueptr, zbufptr, countbufptr, cuezbufptr, sumbufptr); SetPort(tPort); DestInfo := info; info := SourceInfo; SelectSlice(info^.StackInfo^.CurrentSlice); if AutoSelectAll then KillRoi; info := DestInfo; end; {procedure DoProjection} {******************************************************************************} {* This procedure allows the user to set parameters for projection in a large dialog box. *} {******************************************************************************} function ShowProjectDialogBox: boolean; const ProjectOptionsID = 128; SliceIntervalID = 3; InitAngleID = 4; TotalAngleID = 5; AngleIncID = 6; TransparencyLowerID = 7; TransparencyUpperID = 8; OpacityID = 9; DepthCueSurfID = 10; DepthCueIntID = 11; RotationAboutXID = 12; RotationAboutYID = 13; RotationAboutZID = 14; SaveProjectionsID = 15; MinProjSizeID = 16; NearestID = 28; BrightestID = 29; MeanID = 30; var mylog: Dialogptr; {pointer to dialog box} i, item, alldone: integer; SaveInitAngle, SaveTotalAngle, SaveAngleInc: integer; SaveOpacity: integer; SaveAxisOfRotation: AxisType; SaveSaveProjections, SaveCloseSlices, SaveMinProjSize: Boolean; SaveLower, SaveUpper: integer; PercentSurf, PercentInt: integer; interval: extended; begin if ForceToFront <> noErr then begin ShowProjectDialogBox := false; exit(ShowProjectDialogBox); { AE - RMD 1/10/95 } end; InitCursor; Interval := info^.StackInfo^.SliceSpacing; if Interval <= 0.0 then Interval := 1.0; SaveInitAngle := InitAngle; SaveTotalAngle := TotalAngle; SaveAngleInc := AngleInc; SaveOpacity := Opacity; SaveAxisOfRotation := AxisOfRotation; SaveSaveProjections := SaveProjections; SaveMinProjSize := MinProjSize; SaveLower := TransparencyLower; SaveUpper := TransparencyUpper; PercentSurf := 100 - DepthCueSurf; PercentInt := 100 - DepthCueInt; if DensitySlicing then with info^ do begin TransparencyLower := SliceStart; TransparencyUpper := SliceEnd; end; mylog := GetNewDialog(ProjectOptionsID, nil, pointer(-1)); SetDReal(MyLog, SliceIntervalID, Interval, 1); SelectdialogItemText(MyLog, SliceIntervalID, 0, 32767); SetDNum(MyLog, InitAngleID, InitAngle); SetDNum(MyLog, TotalAngleID, TotalAngle); SetDNum(MyLog, AngleIncID, AngleInc); SetDNum(MyLog, TransparencyLowerID, TransparencyLower); SetDNum(MyLog, TransparencyUpperID, TransparencyUpper); SetDNum(MyLog, OpacityID, Opacity); SetDNum(MyLog, DepthCueSurfID, PercentSurf); SetDNum(MyLog, DepthCueIntID, PercentInt); OutlineButton(MyLog, ok, 16); SetDlogItem(MyLog, RotationAboutXID, ord(AxisOfRotation = XAxis)); SetDlogItem(MyLog, RotationAboutYID, ord(AxisOfRotation = YAxis)); SetDlogItem(MyLog, RotationAboutZID, ord(AxisOfRotation = ZAxis)); SetDlogItem(MyLog, NearestID, ord(ProjectionMethod = NearestPoint)); SetDlogItem(MyLog, BrightestID, ord(ProjectionMethod = BrightestPoint)); SetDlogItem(MyLog, MeanID, ord(ProjectionMethod = MeanValue)); SetDlogItem(MyLog, SaveProjectionsID, ord(SaveProjections)); SetDlogItem(MyLog, MinProjSizeID, ord(MinProjSize)); alldone := 0; repeat {if we don't do it this way, ModalDialog throws us into code checking after each keystroke} repeat {meaning you can't type in a 2 digit number} ModalDialog(nil, item); if item = SaveProjectionsID then begin SaveProjections := not SaveProjections; SetDlogItem(MyLog, SaveProjectionsID, ord(SaveProjections)); end; if item = MinProjSizeID then begin MinProjSize := not MinProjSize; SetDlogItem(MyLog, MinProjSizeID, ord(MinProjSize)); end; if (item = RotationAboutXID) or (item = RotationAboutYID) or (item = RotationAboutZID) then begin case item of RotationAboutXID: AxisOfRotation := XAxis; RotationAboutYID: AxisOfRotation := YAxis; RotationAboutZID: AxisOfRotation := ZAxis; end; {case} SetDlogItem(MyLog, RotationAboutXID, ord(AxisOfRotation = XAxis)); SetDlogItem(MyLog, RotationAboutYID, ord(AxisOfRotation = YAxis)); SetDlogItem(MyLog, RotationAboutZID, ord(AxisOfRotation = ZAxis)); end; if (item >= nearestID) and (item <= MeanID) then begin case item of NearestID: ProjectionMethod := NearestPoint; BrightestID: ProjectionMethod := BrightestPoint; MeanID: ProjectionMethod := MeanValue; end; SetDlogItem(MyLog, NearestID, ord(projectionMethod = NearestPoint)); SetDlogItem(MyLog, BrightestID, ord(projectionMethod = BrightestPoint)); SetDlogItem(MyLog, MeanID, ord(projectionMethod = MeanValue)); end; until (item = ok) or (item = cancel); alldone := 1; if (item = ok) then begin {check all the values that could have been entered, don't know which were changed} Interval := GetDReal(MyLog, SliceIntervalID); if (Interval <= 0.0) or (Interval > 1023.0) then begin Interval := info^.StackInfo^.SliceSpacing; SetDReal(MyLog, SliceIntervalID, Interval, 1); beep; alldone := 0; end; {if Interval} InitAngle := GetDNum(MyLog, InitAngleID); if (InitAngle < 0) or (InitAngle > 360) then begin InitAngle := SaveInitAngle; SetDNum(MyLog, InitAngleID, InitAngle); beep; alldone := 0; end; {if InitAngle} TotalAngle := GetDNum(MyLog, TotalAngleID); if (TotalAngle < 0) or (TotalAngle > 360) then begin TotalAngle := SaveTotalAngle; SetDNum(MyLog, TotalAngleID, TotalAngle); beep; alldone := 0; end; {if TotalAngle} AngleInc := GetDNum(MyLog, AngleIncID); if (AngleInc < 0) or (AngleInc > 360) then begin AngleInc := SaveAngleInc; SetDNum(MyLog, AngleIncID, AngleInc); beep; alldone := 0; end; {if AngleInc} TransparencyLower := GetDNum(MyLog, TransparencyLowerID); if (TransparencyLower < 0) or (TransparencyLower > 255) then begin TransparencyLower := SaveLower; SetDNum(MyLog, TransparencyLowerID, TransparencyLower); beep; alldone := 0; end; {if TransparencyLower} TransparencyUpper := GetDNum(MyLog, TransparencyUpperID); if (TransparencyUpper < 0) or (TransparencyUpper > 255) then begin TransparencyUpper := SaveUpper; SetDNum(MyLog, TransparencyUpperID, TransparencyUpper); beep; alldone := 0; end; {if TransparencyUpper} Opacity := GetDNum(MyLog, OpacityID); if (Opacity < 0) or (Opacity > 100) then begin Opacity := SaveOpacity; SetDNum(MyLog, OpacityID, Opacity); beep; alldone := 0; end; {if Opacity} PercentSurf := GetDNum(MyLog, DepthCueSurfID); if (PercentSurf < 0) or (PercentSurf > 100) then begin PercentSurf := 100 - DepthCueSurf; SetDNum(MyLog, DepthCueSurfID, PercentSurf); beep; alldone := 0; end; {if DepthCueSurf} PercentInt := GetDNum(MyLog, DepthCueIntID); if (PercentInt < 0) or (PercentInt > 100) then begin PercentInt := 100 - DepthCueInt; SetDNum(MyLog, DepthCueIntID, PercentInt); beep; alldone := 0; end; {if DepthCueInt} info^.StackInfo^.SliceSpacing := Interval; end; until (alldone = 1); DisposeDialog(mylog); if item = cancel then begin {if Cancel, keep the saved values} InitAngle := SaveInitAngle; TotalAngle := SaveTotalAngle; AngleInc := SaveAngleInc; Opacity := SaveOpacity; AxisOfRotation := SaveAxisOfRotation; SaveProjections := SaveSaveProjections; MinProjSize := SaveMinProjSize; TransparencyLower := SaveLower; TransparencyUpper := SaveUpper; ShowProjectDialogBox := false; end else begin DepthCueSurf := 100 - PercentSurf; DepthCueInt := 100 - PercentInt; ShowProjectDialogBox := true; end; end; procedure Project; begin if ShowProjectDialogBox then DoProjection; end; end.